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Abstract 

We present a program for the numerical evaluation of form factors entering the 
calculation of one-loop amplitudes with up to six external legs. The program is 
written in Fortran95 and performs the reduction to a certain set of basis integrals 
numerically, using a formalism where inverse Gram determinants can be avoided. 
It can be used to calculate one-loop amplitudes with massless internal particles in 
a fast and numerically stable way. 
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Computer: Any computer with a Fortran95 compiler 
Operating system: Linux, Unix 

RAM: RAM used per form factor is insignificant, even for a rank six six-point form 
factor 

Number of processors used: one 

Keywords: NLO Computations, One-Loop Diagrams, Tensor Reduction 
PACS: 12.38.Bx 

Classification: 4.4 Feynman diagrams, 11.1 High Energy Physics Computing 
External routines /libraries: perl 

Nature of problem: Evaluation of one-loop multi-leg tensor integrals occurring in 
the calculation of next-to-leading order corrections to scattering amplitudes in ele- 
mentary particle physics. 

Solution method: Tensor integrals are represented in terms of form factors and a 
set of basic building blocks ( "basis integrals" ) . The reduction to the basis integrals 
is performed numerically, thus avoiding the generation of large algebraic expressions. 

Restrictions: The current version contains basis integrals for massless internal par- 
ticles only. Basis integrals for massive internal particles will be included in a future 
version. 

Running time: Depends on the nature of the problem. A rank 6 six-point form 
factor at a randomly chosen kinematic point takes 0.13 seconds on an Intel Core 2 
Q9450 2.66 GHz processor. 
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LONG WRITE-UP 



1 Introduction 

Collider experiments at the TeV scale, in particular the LHC experiments, 
are expected to shed light on the mechanism of electroweak symmetry break- 
ing and to open up new horizons concerning our understanding of elementary 
particle interactions. In order to achieve these goals, expected signal as well 
as background rates should be well under control, which implies that a multi- 
tude of scattering processes should be known at next-to-leading order (NLO) 
accuracy. 

Over the last years, enormous efforts have been made to calculate NLO cor- 
rections, in QCD as well as in the electroweak sector. For a review see e.g. pQ. 
These calculations in general involve two parts, the treatment of extra real 
emission and the calculation of virtual corrections, i.e. one-loop amplitudes. 
While the calculation of one-loop amplitudes with up to four external particles 
has reached a quite mature state meanwhile, and automated tools have been 
developed already some time ago PIBIHIqTIo] . the calculation of processes with 
five or more external legs required and boosted new developments in various di- 
rections, for recent developments see e.g. [T1IT1I81IMITM 

Initially, NLO calculations have mostly been done on a process-by-process ba- 
sis, but fortunately we are moving towards automation also for multi-particle 
processes, as can be seen from the tools which have been constructed re- 
cently. For the automated calculation of one-loop amplitudes with more than 
four external legs, there are the publicly available programs FormCalc/Loop- 
Tools |4"f6] which recently have been extended to 5-point processes [24125] 
and the program CutTools [26] which is based on a numerical unitarity for- 
malism [TTfl5|ll9ll20f22j . Further, there are the programs BlackHat [27] and 
Rocket [22], relying also on cutting techniques. Concerning the generation of 
subtraction terms for real radiation, automated tools have become publicly 
available recently as well [29.3 0|l3T|32] . Integral libraries for massive [2|25] as 
well as infrared divergent [33] scalar integrals also exist. 

As already mentioned, a public program for the reduction of tensor integrals 
so far is available only for infrared-finite integrals and for up to five external 
legs [4.25J. In this paper we present a program for the numerical reduction of 
tensor integrals with up to six external legs. In the present version, we focus 
mainly on massless QCD applications, i.e. processes with massless internal 
particles. The master integrals are implemented in the code to be valid in 
all kinematic regions. Infrared divergences are regulated dimensionally, i.e the 
loop momenta live in n = 4— 2e dimensions. The output for a specific kinematic 
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point is a set of six numbers representing the real and imaginary parts of the 
coefficients of the Laurent series in e, i.e. the coefficients of the 1/e 2 , 1/e poles 
and the finite part. 

The reduction formalism is valid for massless as well as massive external and 
internal particles. However, the basis integrals for processes involving internal 
massive particles will be implemented in a forthcoming version. We would like 
to emphasize that the program can be used not only for tensor reduction, but 
also to calculate basis integrals, with or without Feynman parameters in the 
numerator, and therefore is also of interest for calculations where the integral 
coefficients have been determined by unitarity cut techniques: golem95 can be 
used library for master integrals. 

The paper is organised as follows. In section 2, we review shortly the theoreti- 
cal background. Section 3 contains a brief summary of the software structure, 
while section 4 contains a detailed description of the individual components of 
the program. The installation instructions are given in section 5, and section 6 
contains the descriptions of three different test runs: the calculation of a form 
factor for a rank five five-point function, the calculation of all form factors 
for rank one six-point functions in one go, and finally the calculation of a full 
amplitude: the helicity amplitudes for light-by-light scattering. We give an 
outlook on future versions in section 7 and explain technical details in appen- 
dices \Kj] and [A~2] The code comes with a number of demonstration programs, 
demonstrating for example the behaviour near a scattering singularity, or the 
relation to LoopTools notation. All these demonstration programs are listed 
in Appendix IA.31 



2 Theoretical background 

The program is an implementation of the formalism developed in Ref. [7j . Here 
we will summarize only its main features relating to the golem95 program, for 
further details we refer to [7]. 



2.1 Form Factors 

A general one-loop tensor integral of rank r can be written as 
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Fig. 1. General iV-point one-loop graph with momentum and propagator labelling. 



where q a = k + r a , and r a is a combination of external momenta. For the 
diagram in Fig.[H r« = J2)=iPj- Our method is defined in n = 4— 2e dimensions 
and thus is applicable to general scattering processes with arbitrary propagator 
masses. Taking integrals of the form ([T]), i.e. with q£ instead of just k^ in 
the numerator, as building blocks has two advantages: first, combinations of 
loop and external momenta appear naturally in Feynman rules, second, it 
allows for a formulation of the tensor reduction which manifestly maintains 
the invariance of the integral under a shift k — > k + tq in the loop momentum. 
Such a shift can be absorbed into a redefinition of the rj, rj —>■ rj — tq. By 
setting oi, . . . , a r — N, and using momentum conservation to set r^ = 0, we 
can always retrieve the commonly used form 

/rl n h Jet 11 h^ r 
_ ' ' ' (0) 

% W 2 (qf -ml + i5)... (q% - m\ + iS) ' K ' 

The Lorentz structure of the integral ([1]) is carried by tensor products of the 
metric g^ u and the difference vectors 

K = r? - rf , (3) 



which are shift invariant. Therefore, tensor integrals are expressible by linear 
combinations of such Lorentz tensors and form factors A ;^ , B^' r _, C N : r , 
defined by 



A, • • ■ A'- . 

Jl jr- 
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+ V \g-g-A) •••A;, A^"' 1 ** C?' r < (S) (4) 

where [■ • denotes the distribution of the r Lorentz indices /ij, and 

momentum labels a« to the vectors A^ l a . in all distinguishable ways. S denotes 
an ordered set of propagator labels, related to the kinematic matrix S, defined 
by 



Sn = (n ~ T if ~ m t ~ ™% ; i,j e{i,...,N} . (5) 

There is a one-to-one correspondence between «Sy and the set S = {1, . . . , N}. 
We recall that standard form factor representations can be simply obtained 
by replacing a,- = iV for all j, together with = 0. This also shows that 
the form factors do not depend on the introduction of the difference vectors 
Afj. The form factors are shift invariant by themselves. Therefore the program 
golem95 can be used without ever introducing difference vectors, if the user 
prefers not to do so. 

Due to the fact that for N > 5, four linearly independent external vectors 
form a basis of Minkowski space, the tensor reduction for N > 6 can be done 
in such a way that only form factors for N < 5 are needed. Therefore, the 
Lorentz structure of (N > 5)-point rank r tensor integrals does not require the 
introduction of additional factors of g^ u as compared to the N = 5 case, only 
additional external vectors appear. We note that for N = 5, one could already 
express the metric by external momenta, but this would introduce inverse 
Gram determinants. In [34], it is shown that all tensor five-point functions 
can be reduced to some basis integrals without generating higher dimensional 
five-point functions. In [7], a formal proof of this fact can be found, as well as a 
reduction method which avoids both inverse Gram determinants and spurious 
higher dimensional five-point functions. A method where inverse Gram deter- 
minants in the reduction from five-point to four-point integrals are absent is 
also presented in Ref. [9]. 

The form factors are linear combinations of reduction coefficients and ba- 
sid~n integrals, where our basis integrals are not necessarily scalar integrals, 
as explained in section 12.31 The reduction coefficients are derived from the 
kinematic matrices S, where we define 



i E*.- B = E 6 *- ( 6 ) 

fees ies 



1 We call them "basis integrals " because they are the endpoints of our reduction, 
although they do not form a basis in the mathematical sense. 
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The quantity B is related to the Gram determinant by 
B detS = (-l) N+1 detG . 



(7) 



The form factors are all given explicitly in [7]. As an example, a rank two 
pentagon integral is represented as 



J 5 ™(a 1 ,a 2 ;5)= £ Af^ A^S) + ^ B 5 >\S) 
h,heS 

B 5 ' 2 (S) = -W b i I 2 +2 (s\{j}) 



+ ^E E feC+^id^^U*}) 

z jes kes\{j} 

and it is the form factors like A^ 2 l2 (S), B 5 ' 2 (S) which are implemented in 
golem95. 



The program golem95 can be used for amplitude calculations in several ways: 
One approach, which aims to avoid tensor integrals of high rank, is to can- 
cel the reducible numerators of an expression before interfacing to golem95 
to calculate the irreducible tensor integrals. However, as all form factors for 
maximal rank (in a renormalisable gauge) are implemented, the expression 
for an amplitude can also be interfaced to golem95 without performing any 
cancellations between numerators and propagators. This has the advantage 
that the algebraic manipulations to do are minimal, and that even the Dirac 
traces, which often appear as coefficients of the form factors, can be done 
numerically. 



2.2 Feynman parameter representations 



The golem95 program uses the fact that tensor integrals are related to Feyn- 
man parameter integrals with Feynman parameters in the numerator. The 
basic object is the set S, containing the labels of the propagators which define 
the integral. A scalar integral, after Feynman parametrisation, can be written 

as 

„ r N N „_ N 

ius)=(-i) N nN-^) /n 

i=i i=i 
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R 2 = -g E -itf 



(9) 



In general, a one-loop iV-point amplitude will contain iV-point integrals as well 
as (N — 1), (N — 2), . . . , (N — M)-point integrals with tree graphs attached to 
some of the external legs of the loop integral. The latter are characterised by 
the omission ("pinch") of some propagators (say ji, . . . ,j m ) of the "maximal" 
one loop iV-point graph, and therefore correspond to a subset of S where 
certain propagator labels are missing, S \ {ji, . . . ,j m }. The program golem95 
is based on this concept of sets characterising the integrals. 

The general relation between tensor integrals and parameter integrals with 
Feynman parameters in the numerator is well known [35,36,37.38] 
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where I^ +2m (ji . . . , j r -2m ; S) is an integral with Feynman parameters in the 
numerator, [r/2] stands for the nearest integer less or equal to r/2 and the 
symbol ®m indicates that m powers of the metric tensor are present. Feynman 
parameter integrals corresponding to diagrams where propagators ji , . . . , j m 
are pinched with respect to the "maximal" topology can be defined as 
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2.3 Basis integrals 



The basis integrals, i.e. the endpoints of our reduction, are 4-point functions in 
6 dimensions i|, which are IR and UV finite, UV divergent 4-point functions 
in n + 4 dimensions, and various 2-point and 3-point functions, some of the 
latter with Feynman parameters in the numerator. This provides us with a very 
convenient separation of IR/UV divergences, as the IR poles are exclusively 
contained in the triangle functions. Explicitly, our reduction basis is given by 
integrals of the type 
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to., ...,*>= -r (s - =) / n * *a - e ,> H ,.|;:!V^ • 
/r o-o = -r (2 - =) / n * «d - £ ,> Hz . s :':_^ n , 
iruu ■ ■ ■ Jr) = r (3 - 1) J n *a-x>) H ..£;V* • 
™ =r (2 - \) /n ^ - X>) ( _, 2 . 5 .:_ i , )2 -„ /2 • 

(12) 

where r max = 3, as well as I£ , I^ +2 , 12 +2 , 12 + with no Feynman parameters 
in the numerator, and two-point functions. 

Note that I^ +2 and J4" 1 " 4 are UV divergent, while J3 can be IR divergent. In 
the code, the integrals are represented as arrays containing the coefficients of 
their Laurent expansion in e = (4 — n)/2. 

Further reduction of these integrals to scalar basis integrals (i.e. integrals 
with no Feynman parameters in the numerator) introduces factors of 1/B, 
i.e. inverse Gram determinants. A particular feature of golem95 is the fact 
that the above integrals are not reduced to scalar basis integrals in cases 
where B becomes small, thus avoiding problems with small inverse deter- 
minants. In these cases, the above integrals are evaluated numerically. As 
B = (— 1) N+1 det(G?)/ det(iS) is a dimensionful quantity, the switch to the nu- 
merical evaluation of the basis integrals is implemented such that the value of 
the dimensionless parameter B is tested, where 

B = B x (largest entry of S) . (13) 



If B > B cut , the reduction is performed, else the program switches to the 
direct numerical evaluation of the integral. The default value is B cut = 0.005. 
A major improvement with respect to the numerical evaluation method used 
in [7] is the following: while in [7] the numerical evaluation of box integrals 
was based on three-dimensional parameter representations, we use a certain 
one-dimensional parameter representation here, obtained after performing two 
integrations analytically, as outlined in Appendix IA.2I In this way one can 
use deterministic integration routines, leading to a fast and precise numerical 
evaluation. This has been done for box integrals with up to three off-shell legs 
and all triangle integrals. The relative error to be achieved in the numerical 
integration has been set to the default value 10 -8 . If this precision has not 
been reached, the program will write a message to the file error.txt. In 
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some cases, calculating in double precision Fortran may not be sufficient. The 
code is designed such that it can be compiled in quadruple precision as well. 

We would like to emphasize that the program also can be used as a library 
for master integrals with massless internal particles. For example, the scalar 
box integrals in n dimensions, with up to four off-shell external legs, can be 
calculated by just calling the form factor A 4 ' . Depending on the kinematics, 
the program will call the appropriate box type automatically. The scalar box 
integrals in n + 2 and n + 4 dimensions are related to the form factors by 
S 4,2 = _ /4 «+2/ 2; c 14 - 4 = J 4 n+4 /4, analogous for N = 3. 



3 Overview of the software structure 

The structure of the golem95 program is the following: There are four main 
directories: 

(1) src: the source files of the program 

(2) demos: some programs for demonstration 

(3) doc: documentation which has been created with robodoc [39] 

(4) test: supplements the demonstration programs, containing files to pro- 
duce form factors with user-defined kinematics. The user can specify the 
rank, numerator, numerical point etc. via a steering file. 



4 Description of the individual software components 

Here we give a short summary of the contents of the individual modules. A 
detailed description of the usage, dependencies and output of each module is 
given at the beginning of each file of the program and can also be read in html 
format by loading the file masterindex.html from the subdirectory doc into 
the browser and following the various links. 

The program is written in Fortran95 and is downwards compatible to For- 
tran90. 

The directory src contains the subdirectories 

• form_factor: contains five modules to compute the form factors for two- 
point to six-point functions: 

f orm_f actor_2p.f 90, f orm_f actor_3p . f 90 , f orm_f actor_4p . f 90, 
f orm_f actor_5p . f 90 , f orm_f actor_6p . f 90. 

• integrals: contains the subdirectories four _point, three_point, two_point. 
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four_point: contains six modules to compute the four-point functions with 
pf 7^ holding for four, three, two, one or none of the external legs: 
function_4plm.f 90, function_4p2m_opp.f 90, f unction_4p2m_adj .f90, 
function_4p3m.f 90, function_4p4m.f 90, generic_f unction_4p . f 90. 

three_point: contains six modules to compute the three-point functions 
with three, two or one external legs off-shell: 

function_3plm. f 90, function_3p2m. f 90, function_3p3m. f 90, 
generic_function_3p.f 90, mod_h0.f90, mod_hf.f90, mod_he.f90. 
two_point: contains one module to compute the two-point functions: 
generic_f unction_2p . f 90. 

• kinematic: contains two modules to compute the matrix S and its inverse 
and to compute the reduction coefficients be 

matrice_s . f 90 , inverse_matrice . f 90. 

The definition of B (see eq. ( 1131) ) is contained in matrice_s . f 90. 

• module: contains auxiliary functions/subroutines and the definition of 
some default parameters: 

The file parametre . f 90 contains the parameters defining the switch be- 
tween the reduction down to scalar basis integrals (which are implemented 
in analytic form) versus the numerical evaluation of integrals (with or with- 
out Feynman parameters in the numerator), as explained in section 12.31 
The default value for B for three-point as well as four-point functions has 
been set to 0.005. The other default parameters for the numerical integra- 
tion are also fixed in parametre . f 90. Further, there is a switch to calculate 
the rational parts of amplitudes only. The default is tot to calculate the 
complete form factors. If tot is replaced by rat, only the rational parts will 
be calculated. 

The auxiliary functions will not all be listed here, we only point to the most 
important ones: 

• Polylogarithms and other special functions are defined in z_log.f90, 
zdilog.f90, kronecker.f 90, constante . f 90. 

• spinor . f 90 contains functions to compute scalar products of four-momenta, 
spinorial products and totally antisymmetric epsilon tensors. 

• The files preci_double . f 90 and preci_quad . f 90 are needed for the switch 
between double precision and quadruple precision. The default is dou- 
ble precision. If quadruple precision should be used, one has to define 
$precision = "quadruple" in the file configure.pl. Note that quadru- 
ple precision is at present only supported by the if ort compiler. 

• The file f orm_f actor_type . f 90 defines a type form-factor such that form 
factors, which are arrays of three complex numbers, can be involved in 
algebraic manipulations. 

■ c ache. f 90 is used to reserve memory in order to store results for three- 
point or four-point functions which already have been computed. 

• numerical: contains two modules for the numerical integration : 
mod_adapt_gauss . f 90 , mod_numeric . f 90. 
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Concerning the numerical integration, the following features should be pointed 
out: 

• The user can change the integration method for the numerical integra- 
tion of the one-dimensional parameter integrals by changing the module 
numerical_evaluation in the file mocLnumeric . f 90 in the directory 
src/numerical. 

• The values for the cuts defining the switch to a one-dimensional numerical 
integration of the basis integrals are given in parametre . f 90 and can be 
changed easily by the user. 

Note: If the user wants to change the default values defined in parametre . - 
f 90, it is not necessary to recompile the library If the desired values are 
defined in the main program, the default values will be overwritten. The 
command use parametre still has to be included in the header of the main 
program. 

• For boxes with 4 off-shell external legs, the expressions for one-dimensional 
numerical integrations are not worked out in this version. Here the program 
will always reduce numerically to scalar basis integrals, irrespective of the 
size of the Gram determinants. 



5 Installation instructions 



The program can be downloaded as a .tar.gz archive from the following URL: 



http : //lappweb . in2p3 . f r/lapth/Golem/golem95_vl . . tar . gz The instal- 



lation instructions given below also can be found in the Readme file coming 
with the code. 



To install the golem95 library, type the following commands: 

./ configure.pl [ — install_path=mypath] [ — compiler=mycompiler] 

make 

make install 



Please note that mypath must be the absolute path of the directory where you 
would like the library to be intalled. If no option for install_path is given, a 
subdirectory of the current directory with the name libgolem will be created 
and the library will be installed in this subdirectory. 



For example, if you want to put the library into the directory 

/home/myname/lib/libgolem and use the compiler g95, then type: 

./ configure.pl — install_path=/home/myname/lib/libgolem — compiler 

make 

make install 
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The directory /home/myname/lib/libgolem will then contain a collection of 
files of type .mod plus a file named libgolem.a which is the golem95 library. 

If no option for the compiler is specified, the installation script will search for 
fortran 95 compilers installed on your system and will take the first matching 
compiler found. 

The program has been tested with the GNU compilers g95 and gfortran, the 
intel compiler ifort, the dec compiler f95, the NAG compiler f95 and the port- 
land compiler pgf95. 



6 Test run description 

The program comes with several demonstration programs located in the sub- 
directory demos. A list of all options contained in the demos directory is given 
in Appendix IA.31 We will describe some selected examples in the following. 

6. 1 Rank five five-point form factor 

As an example for a test run, we first describe the calculation of a form 
factor for a rank five 5-point integral, A^f - 5 . We choose ji = i, i.e. z\ ... z§ in 
the numerator. Further, we choose the following numerical point (in terms of 
entries of the kinematic matrix S containing the invariants): 
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These values are already implemented in the file demo_5point . f 90 in the 
subdirectory demos. All the user has to do is the following: 

• go to the subdirectory demos 

• type "perl configure.pl". The shell will prompt for the choice of the demo to 
be run: 

Choose which demo program you want to run: 

1) three-point functions 

2) four-point functions 



13 



3) five-point functions 

4) six-point functions 

5) 4-photon helicity amplitudes 

6) numerical stability demo: detG — >0 

7) numerical stability demo: det 5 — > 

8) Golem <-> LoopTools conventions 

• Choosing option 3 will produce the following output: 
you have chosen option 3: five-point functions 
The Makefile has been created 

Please run: 
make 

. / comp . exe 

• Running "make" will produce the executable comp . exe where the demo* . f 90 
files matching the choice above will be compiled automatically. Running 
comp . exe will prompt for the rank of the form factor to be calculated: 
Choose what the program should compute: 

0) form factor for five-point function, rank 

1) form factor for five-point function, rank 3 (zl*z2*z4) 

2) form factor for five-point function, rank 5 (zl*z2*z3*z4*z5) 

3) form factor for diagram with propagator 3 pinched, rank 

4) form factor for diagram with propagators 1 and 4 pinched, rankO 

• Choosing option 2 will produce the result which will be written to the file 
test5point.txt and looks as follows: 



The kinematics is: 

Pi + P2 + Pz + Pi + Ps = 
S(l,3) = (p 2 +p 3 ) 2 = -3. 
S(2,4) = (p 3 +p 4 ) 2 = 6. 
5(2,5) = {pi+P2? = 15. 
5(3,5) = ( Pi +p,) 2 = 2. 
5(l,4) = (j9 1 +j9 5 ) 2 = -4. 
5(1,2) =p 2 2 = 0. 
5(2,3) =p 2 3 = 0. 
5(3,4) =pl = 0. 
5(4,5) =p 2 5 = 0. 
5(1,5) =pl = 0. 

A factor T(l + e)r(l - e) 2 /T(l - 2e) (4tt /j 2 ) 6 is factored out from the result. 
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result= 1/e 2 * (O.OOOOOOOOOOE+00 + I* O.OOOOOOOOOOE+00) 

+ 1/e * (O.OOOOOOOOOOE+00 + I* O.OOOOOOOOOOE+00) 
+ (-.8615520644E-04 + I* 0.1230709464E-03) 

CPU time= 7.999000000000001E-003 

We recall that we use the integral measure as in eq. ([TJ. The factor r? = 
T(l + e)r(l — e) 2 /r(l — 2e) (47r fi 2 ) e has been extracted from the integrals to 
comply with the conventions of ref. [7j and the MS subtraction scheme. Note 
that it may be advantageous to call golem95 with rescaled, dimensionless 
invariants (e.g. Sij/fi 2 ), for example in cases where most of the invariants have 
very small numerical values in all kinematic regions. 

6.2 Calculating all possible numerators at once 

In the previous example, the numerical point (and the type of numerator for 
tensor integrals) has been fixed in the demo programs. If the user would like 
to give the the numerical point and the Feynman parameters in the numerator 
as an input, he can use the file param. input in the subdirectory test. This 
setup also allows to calculate all possible numerators for a certain rank in one 
go. A typical example looks as follows: 

Assume we would like to calculate the form factors for rank one six-point 
functions for all possible numerators Zj,j — 1 . . . 6, for the following numerical 
point (pi = (Ei,Xi,yi,Zi)): 

j9i = (0.5,0.,0.,0.5) 
p 2 = (0.5,0.,0.,-0.5) 

p 3 = (-0.19178191, -0.12741180, -0.08262477, -0.11713105) 
p 4 = (-0.33662712, 0.06648281, 0.31893785, 0.08471424) 
p 5 = (-0.21604814, 0.20363139, -0.04415762, -0.05710657) 

5 

i=i 

To use these momenta, go to the subdirectory test and editf^lthe file momenta. - 
dat, writing each component of the above momenta into a single line. To cal- 
culate an iV-point function, the program will use the first N momenta found 
in momenta.dat (respectively the momenta file specified in param. input). 

2 Alternatively, random momenta can be generated using the program 
mom_rambo.f, adapted from [40J and also contained in the subdirectory test. 
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For iV = 5 and N — 6, it is important that momentum conservation, i.e. 
Y,iLiPi — 0, is fulfilled because momentum conservation has been assumed to 
hold true in the reduction. 

To generate results, the user only has to do the following: 

• edit the file param . input to choose the number of legs, rank and numera- 
tor. If only a particular numerator should be calculated, give the labels of 
Feynman parameters, else put all into the numerator field, 

• type perl maketest.pl. 

The program will automatically compile the corresponding functions and run 
the executable. The following output will be produced: 

(1) a separate output file called N[nb of legs] [rank] _[pt]. out for each individual 
numerator ([pt] denotes the "label" of a particular numerical point, which 
can be chosen by the user to distinguish results for different numerical 
points) 

(2) a file called N[nb of legs] [rank] _[pt]. numbers, where all form factors that 
have been calculated for a particular rank and number of legs and numer- 
ical point are appended. For example, if the option all has been chosen to 
calculate all possible combinations of Feynman parameters in the numer- 
ator, this file will contain the results for all these numerators. The format 
is such that it can be read by Mathematica, to allow direct comparisons to 
results obtained from algebraic programs. If the result is P 2 /e 2 + P 1 /e + P Q 
for a rank r iV-point form factor of type A, the output will be a list 
a N \j u ...J r ] = {Ke^Im^KeiP^ImlP&KelPolImiPo]}. 

For example, for rank one six-point functions at the numerical point given 
above ("ptl"), having chosen all in param. input to calculate all six possible 
numerators, the program produces seven output files: N6ranklzi_ptl .out for 
i — 1 . . . 6 and the file N6rankl_ptl .numbers. While the files N6ranklzi_ptl . out 
contain, in addition to the result for the particular numerator, also the kine- 
matic point and CPU time information, the file N6rankl_ptl .numbers just 
lists the results. Note that for N > 5, individual form factors are not uniquely 
defined because the metric tensor can be expressed by external momenta, 
such that individual terms can be shifted between form factors of type A, B 
or C. 

6.3 Calculation of the ^-photon helicity amplitudes 

In order to show how golem95 can be embedded into the calculation of full one- 
loop amplitudes, we give here the calculation of the light-by-light scattering 
amplitude in massless QED as a pedagogical example. 
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This amplitude is denned by six Feynman diagrams where the four photons 
are attached to a closed fermion loop in all possible ways. Diagrams which 
differ by the charge flow only lead to the same value, which leaves us with 
three different topologies defined by the photon orderings 1243, 1234 and 1324, 
respectively. Each diagram is IR finite and UV divergent. The UV divergence 
only cancels in the sum of the diagrams. The results for the three independent 

helicity amplitudes + + ++, + + H — , + H are well known, see for example 

[4T|4*2] . For completeness we list the analytic formulae, omitting the irrelevant 
phases 



„4 ++++ = 8 , ^+++- = -8, 

,4.. r t-u, /t\ t 2 + u 2 / ft\ 2 n\l , i 

^- = -8[l + — log(-) + — (log(-) . (15) 



The analytic expressions in terms of form factors and Mandelstam invariants 
which are given in the demo program demo_4photon. f 90 were obtained as 
follows. After working out the trace of gamma matrices one finds for each 
graph a polynomial in scalar products of polarization vectors, Ej, external 
momenta pj and the D = 4 — 2e dimensional loop momentum k. All reducible 
scalar products, i.e. those which can be written in terms of inverse propagators, 
were cancelled directly. The remaining expressions, containing only irreducible 
scalar products, are proportional to tensor integrals which are transformed 
to form factors using eq. (J4]). Each form factor now has scalar coefficients 
containing polarisation vectors and external momenta, i.e. £j • Ej, Ei • pj, s — 
2pi • P21 t — %P2 ' P3 and u = 2pi ■ p%, where we defined all external momenta 
as incoming. Using spinor helicity methods one can map these coefficients 
to polynomials in the Mandelstam variables s, t and u = —s — t. Choosing 
reference momenta p%, p\, p±, p 3 for the polarization vectors E\, 62, £3, £4 
respectively, one easily can show the following relations [43J relevant for the 
+ + ++ amplitude 



2s 2 

et- £ t = ~-r £ t ■ P3 £ t ■ Ps , et-et = 7 £ i ' £ t ' Pi > 
tu t 

2 2,s 
£+•£+ = -£+ • p 3 et ■ pi , = — £+ ■ p A ej ■ p 2 , 

2s 2s 

+ + + + + + -4- + 

£2 - £ I = — £ 2 -Pssj -P2 , ej -ej = -—£3 -Piel-px, 
e + -v e+-v e + -v e + ■ v -KM. 

e x Pse 2 p 3 e 3 Pl e, Pi-[ Ys ) (12)(34) " (16) 

The phase factor in the last line is irrelevant for observables and thus can be 
dropped. All coefficients of the form factors of the + + ++ amplitude are now 
rational polynomials in s, t and u = — s — t. For the + + H — amplitude one 
needs instead of eq. ([T5]) 
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ej -e A = -e{ ■ p 4 el ■ p 1 , ej ■ e 4 = -ej • p 4 ■ p 2 , 

4 • £ 4 = o 

£l - P3 ' 2 ' Pl£ * - Pl -{zs) (12) [41] (13) (17) 

and for the + H amplitude 



2 , 2 

£l ■ ^3 = ~ £ 1 ' 4 • Pi > £ t ■ £ 4 = 7 £ 1 • P4 £4 • Pi > 
4 • £ 3 = • P3 £3 • P2 , £ 3 ■ £ 4 = ~^ £ 3 ' ^ £ * ' P 1 



tw\ 2 [21] (34) 



e 1 -Pse 2 .p s s 3 . Pl s A . Pl =^- J (12)[34] . 

These relations define the different coefficients of the form factors present in 
the file demo_4photon.f 90 which evaluates the four photon amplitude^ 3 "! The 
form factors for each momentum ordering have to be evaluated only once. 
We have compared our numerical result with the well-known results for these 
amplitudes and find perfect agreement. 

The program demo_4photon . f 90 can be used as a guideline how to express any 
amplitude with massless loops in terms of form factors and scalar coefficients 
before evaluating it with golem95. 



7 Conclusions and Outlook 



We have presented the Fortran 95 program golem95 for the numerical eval- 
uation of tensor integrals up to rank six six-point functions. The program is 
based on a form factor representation of tensor integrals and performs the 
reduction to a certain set of basis integrals numerically. The basis integrals 
are implemented in analytic form. If during the reduction process an inverse 
determinant becomes small, the program switches to a numerical evaluation of 
the (tensor-) integral without further reduction, thus avoiding small denomina- 
tors. The numerical evaluation is based on one-dimensional parameter integral 
representations for most of the basis integrals, allowing for a fast and precise 
numerical integration. 



3 We note that the three helicity amplitudes can be evaluated in a much simpler 
way. By applying spinor helicity methods at an earlier stage, one can achieve a 
representation without any tensor four-point function. The given representation 
should only illustrate a generic form factor representation of an amplitude. 
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The results are given as a set of three complex numbers representing the 
coefficients of the Laurent series in the dimensional regularisation parameter 
e, i.e. the coefficients of the 1/e 2 , 1/e poles and the finite part. 

The program can also be used as a library for master integrals (including 
infrared divergent ones) , as the form factors with no Feynman parameter labels 
directly correspond to scalar integrals. In the current version, master integrals 
with massive internal particles are not implemented yet. They will be available 
in a forthcoming version. There is no restriction on the number of massive 
external legs. 

A future version will also combine the golem95 code for the form factor eval- 
uation with a code for the generation of amplitudes, thus moving towards a 
full automatisation of the calculation of one-loop amplitudes. 
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A Appendices 

A.l Landau singularities 

Besides the spurious appearence of powers of inverse Gram determinants 
caused by the decomposition on scalar integrals in the reduction process, 
which can be avoided e.g. with the method advocated here and in [7], an- 
other source of problems in the numerical evaluation of scattering amplitudes 
may be caused by the occurence of actual kinematic singularities, the so-called 
Landau singularities [H]. The latter may appear in some diagrams contribut- 
ing to the considered amplitude whenever the determinant of the kinematic 
matrix S associated with these diagrams - or with reduced diagrams obtained 
by one or several pinches - vanishes. Typical cases of such singularities are 
threshold singularities in loop calculations with internal and external masses, 
collinear and infra-red singularities with massless internal and external lines. 
Another type is the one of scattering singularities PPSS], for which (det G) — > 
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and det S becomes proportional to (det G) 2 , such that both vanish simultane- 
ously. The occurence of these particular cases of vanishing Gram determinants 
should not be confused with the spurious ones. Individual diagrams lead to 
infinities at such kinematic configurations, where a mass singularity and a 
scattering singularity coincide. 

In addition, it should be noted that for a given diagram, the reduction al- 
gorithm breaks down3 at such a scattering singularity, as inferred from the 
relation B oc det Gj det S — > oo. As one combines the diagrams to scattering 
amplitudes, gauge cancellations may occur analytically, which in general re- 
duces the degree of singularity as compared to individual diagrams, or even 
make the singularity bounded, as e.g. observed in the 6 photon amplitudes^! 
[T1H7] . On the other hand, the numerical combination of the singularities from 
separate diagrams is expected to be problematic, and leads to instabilities even 
in cases of expected finiteness. Note that this problem is known to all methods 
based on the reduction of diagrams, so it is is not specific to our reduction for- 
malism. We note that the problem of large numerical cancellations is inherent 
to any method based on the reduction to scalar master integrals like 1^,12, 
as the latter may become linearly dependent near such singularities. 

Depending on the inclusiveness of the observable to be calculated, and the 
degree of the singularity, possible cures could be to resort to mutliple precision 
in some vicinity of the kinematic singularities, and/or place a hole in the phase 
space around the singularity together with a smooth interpolation over it. 
Certainly, in specific cases, when the observable would be controlled by the 
singularity, such methods would be inadequate. 



A. 2 One- dimensional integral representations 



In this appendix we will derive representations of IR finite box-and triangle 
integrals as one-dimensional Feynman parameter integrals. These representa- 
tions have the advantage that they can be integrated numerically in a very 
fast and precise way using deterministic numerical integration routines. This 
approach is similar to the one in jlS] where one parameter integration has 
been carried out analytically. The program switches to this numerical evalu- 

4 In the example of double parton scattering for 2(massless) — > 2 massive legs with 
no internal mass, it can be checked that for the four leg diagram with two opposite 
masses, the equations determining B and 64 in the notations of ref. [7] have no 
solution, because (5v).H.(5v) = 0, hence the equation for B becomes x B = 1 
which has obviously no solution. For a discussion of the reduction in exceptional 
kinematic configurations see also [46J. 

5 Singularity means non-analyticity; the latter can be either infinite - integrable or 
not - or bounded. 
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ation if B becomes smaller than a value defined in module/parametre . f 90 
(the default is 0.005). 



A. 2.1 Four-point functions 

Our starting point are the higher dimensional four-point functions 12 +2 , 12 +4 
given by 



te , . . . , *> = r (3 - 1) ) n ^ «(i-E «> (Ztt^tz^ 

™ = r (2 - 1) / n *<i - g ,) R -^_^ 75 

The reduction of these integrals to integrals with no Feynman parameters 
in the numerator introduces inverse Gram determinants. Therefore it can be 
advantageous to evaluate these integrals without further reduction. To do so, 
we proceed as follows: 

First, to get rid of the 5 distribution, we make the following change of variables: 



Zi = w (1 — x) , z 2 = wxy z , z 3 = wxy (1 — z) , z± = w x (1 — y) (A.l) 

Now, instead of computing directly the three-dimensional integrals numeri- 
cally as proposed in [7], we perform analytically the integration over x and 
y and integrate numerically over the leftover variable z, using an adaptive 
Gauss-Kronrod method 



For the cases treated in the golem95 library (no internal masses), the x and 
y integration for the six- and eight-dimensional four-point functions can be 
computed using two basis integrals: 



dx f =J(n,A,B) , (A.2) 
A + B x 
o 

i 

J dxx n \n(A + Bx)=K(n,A 1 B) (A.3) 
o 

which obey to the following relations: 



J(n, A,B) = ^--^J(n- 1, A, B) (A.4) 
n B B 
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mA , B) J^lzM^} (A . 5) 

K { n, a, B ) = " 1|i ^';;r ( "" u,fl| 

' (A.6) 



(n + 1) 2 

mAB)J A + B > HA + B)-AHA)-B 

B 

Here we assume that B ^ 0; if B = the integrations are trivial. When the 
two first integrations have been done, we are left with the z integration. 

To explain how we proceed, we treat the case of the six-dimensional three- 
mass four-point function as an example. After integration over x and y, we 
are left with the following structure: 



-h \n(h) + e ln(e) h \n(h) - c ln(c) 



fg 



fd 



(A.8) 



with 



c = zSyi + (1 - z) «Si 3 
/ = z (S 24 - S 12 ) + (l-z) (S u - S 13 ) 
g = z(l - z) S23 - zS 2 a - (1 - z)S M 
d = z(l - z) S23 - zS 12 - (1 - z)S 13 
e = zS 24: + (1 - z) (S 3 4 

h = z{l-z)S 23 (A.9) 

where Sij are the iS-matrix elements, they must be understood as Sij + i5 . 

The first thing to note is that / has no poles. All six-dimensional four-point 
functions are infrared finite, and the UV pole of the eight-dimensional four- 
point functions is contained in the overall T-function. Indeed, it is easy to see 
that: g = h — e, d=h — c and / = e — c, so when g — > or d — > or / — > 0, 
the numerator of I goes to zero. To compute the z integral numerically, we 
use a contour deformation: we complexify the z variable 

z = u — ieg(u) (A. 10) 

i.e. we have to compute the following integrals: 
1 1 

J dzf(z) = J duC f(u-ieg(u)) (A.ll) 
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where C is the jacobian of the transformation : C = 1 — i e dg/du and e = ±1. 
The function g has the following properties: g(0) = g(l) = and g(u) > 
for u G [0,1]. For practical applications, we took g{u) — u(l — u). As the 
numerator of / contains some logarithms, some care has to be taken in order 
to avoid a clash between the cut of the logarithm and the contour. To analyse 
that, let us consider the following example: 

E=jdz Ha \ bZ + iS \ X) (A.12) 
J c + dz + is2A 



with a, b, c and d E R, Si, s 2 = ±1 and A > 0. Making the change of variable 
(lAlOl) , we get: 

E= f duC Ha + bu + i (s 1 X-beg(u))) 
J c + du + i (s 2 A — de g(u)) 



By choosing e — — S\ sign(fe), the imaginary part of the argument of the log- 
arithm will be constant and have the sign of s\. This choice of e defines the 
contour but the important point is that by varying u (walking on the contour) 
the cut of the logarithm is never crossed. The pole is located at: 

*o = -~-^A (A.14) 



Using the Cauchy theorem, we arrive at the following relation: 



In (a + bu + i (si X — beg{u))) 



J dzf{z)= J duC- 



c + du + i (s 2 A — de g(u)) 



2i7ri2ee(l + ^)e(-^)«^ (A.15) 



o 

d) 6 i~d) 

where R is the residue of f(z) at z = z . This is the way we proceed to 
compute numerically the two terms of eq. (IA.8I) . We compute the two terms 
separately despite the fact that each term has a pole when / — > while the 
sum does not, because they contain two kinds of logarithms (ln(e) and hi(c)), 
and there is no reason that the choice for e for one term prevents the contour 
from crossing the cut for the other term. 

For the case where there are Feynman parameters in the numerator, every- 
thing works like the preceding example: we always split the integrand of the z 
integration into two pieces (each piece having more terms than the scalar case) 
by separating the two kinds of logarithms. For the other types of four-point 
functions, we proceed in an analogous way. 
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A. 2. 2 Three-mass three-point functions 



In the case of the three-point functions with three off-shell legs, after making 
a change of variables of type flA.ip . we are left with two-dimensional integrals. 
One parameter is integrated out analytically using (1A.2[) . the remaining inte- 
gral is computed numerically, using the same techniques as for the four-point 
functions. 



A. 2. 3 Two-mass three-point functions 

The two mass three-point functions are written in terms of functions Hi [TJ, 
which are defined such that in the numerically problematic case where X —>■ Y , 
their evaluation is numerically stable. The functions Hq, Hi, Hi, H% and H± 
are given by: 



H (X,a) = ^ (A.16) 
H 1 (X,Y,a)= x _ y (A. 17) 

•ya 1 -yl+a 

H 3 (X,Y, a) = — + — w - Jf 

2 Y 2+a -X 2+a 
+ (l + a)(2 + a) (F-X) 3 ( ' 19) 

Y a 3 Vl+a g "p2+a 

HJX, Y, a) = 1 — + — — 

V ' ' ; F-X 1 + a (F-X) 2 (1 + a) (2 + a) (F-X) 3 

g -p3+a _ X^+a 

+ (1 + a) (2 + a) (3 + a) (F - X) 4 (A ' 20) 
X = -X - i 5 

For each function H^X, Y, e), one can define 

Hi(X, Y, e) =eH Ei {X, Y) + ^ tf Fi (X, F) , (A.21) 



and one can show that 



gg ,(X, y)= /^-> z * (A.22) 
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™™-l« ^? + T$? (A - 23) 

v ' 

From this definition, it is easy to show that 



X-Y \n 



±--YH En _ 1 (X,Y)) 



(A.24) 



The equations (IA.22I) and (IA.23I) are used to compute numerically the func- 
tions H En and H Fn . 



A. 3 Contents of the demonstration programs 



The demo programs calculate the following examples, listed also in the file 
DemoContents in the subdirectory demos: 

(1) three-point functions 

(2) four-point functions 

(3) five-point functions 

(4) six-point functions 

(5) calculation of 4-photon helicity amplitudes 

(6) numerical stability demo: det G — > 

(7) numerical stability demo: det S —>■ 

(8) Golem <-> LoopTools conventions 

The items above contain the following options: 



• Three-point functions: 

(1) one off-shell leg 

(2) two off-shell legs 

(3) three off-shell legs 

For each of the three options above, one can choose to calculate: 

(a) scalar three-point function in n dimensions 

(b) three-point function in n dimensions with one Feynman parameter 
(zi) in the numerator 

(c) three-point function in n dimensions with two Feynman parameters 

{zi 22) 

(d) three-point function in n dimensions with three Feynman parameters 

(*?*s) 

(e) scalar three-point function in n+2 dimensions 

(f ) three-point function in n+2 dimensions with one Feynman parameter 

• Four-point functions: 
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(1) no off-shell leg 

(2) one off- shell leg 

(3) two opposite off-shell legs 

(4) two adjacent off-shell legs 

(5) three off-shell legs 

(6) four off- shell legs 

For each of the five options above, one can choose to calculate: 

(a) scalar four-point function in n dimensions 

(b) four-point function in n dimensions with one Feynman parameter (zi) 

(c) four-point function in n dimensions with two Feynman parameters 

(zi z A ) 

(d) four-point function in n dimensions with three Feynman parameters 

(e) four-point function in n dimensions with four Feynman parameters 

\Z\ Z 2 Z3 Z4) 

(f) scalar four-point function in n+2 dimensions 

(g) four-point function in n+2 dimensions with two Feynman parameters 
Oi z 2 ) 

(h) scalar four-point function in n+4 dimensions 

• Five-point functions: 

(1) form factor for five-point function, rank 

(2) form factor for five-point function, rank 3 (z\ z 2 Z4 in numerator) 

(3) form factor for five-point function, rank 5 {z\ z 2 z% z^ z§ in numerator) 

(4) form factor for a pinched 5-point diagram (propagator 3 missing), rank 

(5) form factor for a doubly pinched 5-point diagram (propagators 1 and 4 
missing), rank 

• Six-point functions: 

(1) form factor for six-point function, rank 

(2) form factor for six-point function, rank 4 (zf z 2 z% in numerator) 

(3) form factor A5 for pinched diagram, propagator 3 missing, rank 

(4) form factor for double pinched diagram, propagators 2,5 missing, rank 

(5) form factor for triple pinched diagram, propagators 2,4,6 missing, rank 

• Calculation of 4-photon helicity amplitudes: 

the purpose of this example is to demonstrate how to use golem95 for the 
calculation of full amplitudes. It calculates three different helicity configu- 
rations of the on-shell 4-photon amplitude for a certain kinematic point. 

• Numerical stability demo: det G — » 0: 

calculates a rank three four-point function (in 6 dimensions) in a region 
where \B\ = det Gj det S becomes small, i.e. where a representation based 
on the reduction to scalar integrals would fail. The Feynman parameters 
in the numerator are Z\ z\. The example follows closely the one described 
in section 7.2 of [7J and is also described in the golem95 manuscript: The 
program makes 30 iterations where B = — det Gj det S becomes smaller 
in each iteration. The results for real and imaginary parts of lf(ziz 2 ) are 
written to the file demo_detG.dat as a function of x, where \B\ x 2 for small 
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x. The files plotDetG_Re . gp and plotDetG_Im.gp can be used to plot the 
result with gnuplot by load 'plotDetG_Re/Im.gp' . One can see from the 
plots that The file demo_detG . txt contains the details of the kinematics for 
each iteration. 

• Numerical stability demo: det S — > 0: 

tests the rank 5 five-point tensor coefficient A 5,5 (l, 1, 1, 1, 1) with respect 
to its behaviour when a sub-determinant det S ~ (detG) 2 — > 0. The re- 
sults for real and imaginary parts of the e° part of A 5 ' 5 are written to the 
file demo_a55_dets_sing.dat as a function of the transverse momentum of 
particle 5 and can be plotted with gnuplot by load 'plot_demo_A55.gp'. 

• Relation between Golem output and Loop Tools format: 

produces Golem output for four-point functions up to rank four and gives 
the relation to LoopTools conventions. If LoopTools is linked, the lines con- 
taining the call of LoopTools functions can be uncommented to produce 
LoopTools output in parallel. 
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